A Class of Algorithms for Recovery of Continuous Relaxation Spectrum from Stress Relaxation Test Data Using Orthonormal Functions

The viscoelastic relaxation spectrum provides deep insights into the complex behavior of polymers. The spectrum is not directly measurable and must be recovered from oscillatory shear or relaxation stress data. The paper deals with the problem of recovery of the relaxation spectrum of linear viscoelastic materials from discrete-time noise-corrupted measurements of relaxation modulus obtained in the stress relaxation test. A class of robust algorithms of approximation of the continuous spectrum of relaxation frequencies by finite series of orthonormal functions is proposed. A quadratic identification index, which refers to the measured relaxation modulus, is adopted. Since the problem of relaxation spectrum identification is an ill-posed inverse problem, Tikhonov regularization combined with generalized cross-validation is used to guarantee the stability of the scheme. It is proved that the accuracy of the spectrum approximation depends both on measurement noises and the regularization parameter and on the proper selection of the basis functions. The series expansions using the Laguerre, Legendre, Hermite and Chebyshev functions were studied in this paper as examples. The numerical realization of the scheme by the singular value decomposition technique is discussed and the resulting computer algorithm is outlined. Numerical calculations on model data and relaxation spectrum of polydisperse polymer are presented. Analytical analysis and numerical studies proved that by choosing an appropriate model through selection of orthonormal basis functions from the proposed class of models and using a developed algorithm of least-square regularized identification, it is possible to determine the relaxation spectrum model for a wide class of viscoelastic materials. The model is smoothed and robust on measurement noises; small model approximation errors are obtained. The identification scheme can be easily implemented in available computing environments.


Introduction
Viscoelasticity denotes the joint property of elasticity and viscosity and, hence, describes materials with both fluid and solid properties at the same time. Viscoelastic relaxation or retardation spectra are commonly used to describe, analyze, compare and improve the mechanical properties of polymers [1][2][3][4][5]. The spectra are vital for constitutive models and for the insight into the properties of a viscoelastic material, since, from the relaxation or retardation spectrum, other material functions used to describe rheological properties of various polymers can be uniquely determined [6][7][8][9]. However, the spectra are not directly accessible by measurement. The relaxation and retardation spectra can be recovered from oscillatory shear data and from the time measurements of the relaxation modulus or creep compliance obtained in standard stress relaxation or retardation experiments [1, 2,7]. Many different methods have been proposed during the last five decades for relaxation spectrum computation using data from dynamic modulus tests. Baumgaertel and Winter [10] applied a nonlinear regression for identification of discrete relaxation and retardation time spectra based on dynamic data, in which the number of relaxation times adjusts during the iterative calculations to avoid ill-posedness and to improve the model fit; regularization is not applied here. Honerkamp and Weese [11,12], for relaxation spectrum identification, combined nonlinear regression with Tikhonov regularization and proposed a specific viscoelastic model described by the two-mode log-normal function. Malkin [13] approximated a continuous relaxation spectrum using three constants: the maximum relaxation time, slope in the logarithmic scale and form factor. Malkin et al. [14] derived a method of continuous relaxation spectrum calculations using the Mellin integral transform. An approach proposed by Stadler and Bailly [15] is based on the relaxation spectrum approximation by a piecewise cubic Hermite spline. In turn, Davies and Goulding [16] approximated the relaxation spectrum by a sum of scaling kernel functions located at appropriately chosen points. The algorithm for the relaxation time spectrum approximation by power series was developed by Cho [17]. Anderssen et al. [1] proposed a derivative-based algorithm for continuous spectrum recovery, being also appropriate for the experimental situation where the oscillatory shear data are only available for a finite range of frequencies. These works, but also many others, using different models, approaches, algorithms and computational techniques, have opened new directions of research on discrete and continuous relaxation spectra identification based on dynamic moduli data, which are still being conducted [2,[18][19][20][21].
However, a classical manner of studying viscoelasticity is also through a two-phase stress relaxation test, where time-dependent shear stress is studied for the step increase in strain [4,6,7]. There are only a few papers, e.g., [22][23][24][25][26][27][28], that deal with the spectrum determination from time measurements of the relaxation modulus; additionally, only some of them are addressed to polymers. Therefore, the computationally efficient algorithms to determine the relaxation spectrum applied to time measurements of the relaxation modulus are still desirable. The objective of the present paper was to develop a class of models and an identification algorithm for the continuous relaxation spectrum determination based on discrete-time measurements of the relaxation modulus, which, taking into account the ill-posedness of the original problem of the spectrum recovery, will provide: (a) good approximation of the relaxation spectrum and modulus; (b) smoothness of the spectrum fluctuations, even for noise-corrupted measurements; (c) noise robustness; (d) applicability to a wide range of viscoelastic materials due to the choice of respective model from the considered set of models; (e) ease of implementation of the models and identification algorithm in available computing packages. Thus, the goal of this work was the synthesis of the respective models and general identification scheme, and the analysis of their properties. Approximation errors, convergence, noise robustness, smoothness and the applicability ranges were studied analytically. Further, the numerical verification of the models and algorithm for exemplary theoretical relaxation spectrum, and their applicability to spectrum of real material, polydisperse polymer, was the purpose of this work.
The approach proposed is based on the approximation of the spectrum by a finite linear combination of the basis orthonormal functions. A quadratic identification index, related to the data of the relaxation modulus, is adopted as a measure of the model quality. As a result, the primary infinite dimensional dynamic inverse problem of the continuous relaxation spectrum identification is reduced to the static linear-quadratic programming task. Next, Tikhonov regularization is used to guarantee the well-posed solution. Thus, the approach proposed integrates the technique of an expansion of a function into a series in an orthogonal basis with the least-squares regularized identification [29].
It is demonstrated that due to the choice of appropriate special functions as the basis functions for the unknown relaxation spectrum model, the components in the relaxation modulus model are given by compact analytical or recursive formulas. The technique of expanding an unknown viscoelastic function into a series of orthogonal functions or polynomials has already been used to describe various rheological models of polymers, especially in the time domain. For example, Aleksandrov et al. [30] applied Laguerre polynomials to describe experimentally obtained polyethylene deformation in the creep under diffusion in a liquid environment. Cao et al. [31] used orthogonal expansion based on shifted Legendre polynomials to solve a fractional-order viscoelastic model of polymethyl methacrylate. Abbaszadeh and Dehghan [32] employed a new class of basis based upon the Legendre polynomials to solve a two-dimensional viscoelastic equation. Kim et al. [33] used Chebyshev polynomials for direct conversion of creep data to dynamic moduli.
The idea of using a series of orthogonal functions has also been used to approximate the relaxation spectrum. Lee et al. [34] used the Chebyshev polynomials of the first kind to approximate dynamic moduli data. Stankiewicz [27,28,35] applied orthogonal functions for relaxation spectrum recovery from the stress relaxation data, but these articles use a different definition of the relaxation modulus, according to which the modulus is directly given by the Laplace integral of the spectrum. In this paper, it is shown, for the dominant literature definition of the relaxation spectrum, that the application of the concept of expanding the unknown relaxation spectrum into a series of orthonormal basis functions combined with the least-squares regularized identification allows one to determine the smoothed model of the relaxation spectrum, robust on the measurement noises, with small approximation errors of the relaxation spectrum and modulus. The selection of appropriate orthonormal basis functions, the selection of their time-scale factor and the determination of the optimal regularization parameter using standard generalized cross-validation technique enables the application of the proposed approach to a wide class of viscoelastic materials.

Relaxation Spectrum
The uniaxial, nonaging and isothermal stress-strain equation for a linear viscoelastic material can be represented by a Boltzmann superposition integral [7]: where σ(t) and ε(t) denote the stress and strain at the time t and G(t) is the linear relaxation modulus. Modulus G(t) is given by [1,7,36,37]: or equivalently by where H(τ) and H(v) characterize the distributions of relaxation times τ and relaxation frequencies v, respectively. The continuous relaxation spectra H(τ) and H(v), related by H(v) = H 1 v , are generalizations of the discrete Maxwell spectrum [1,7] to a continuous function of the relaxation times τ and frequencies v. Although other definitions of the relaxation spectrum are used in the literature, for example, in [5,13,28,38], the definition introduced by Equation (2) dominates. The main symbols are summarized in Nomenclature, Appendix C.
The problem of relaxation spectrum determination is the practical problem of reconstructing the solution of the Fredholm integral equation of the first kind (2) or (3) from discrete-time measured data. Time measurements of the relaxation modulus data are considered in this paper. This problem is known to be severely Hadamard ill-posed [39,40].
In particular, small changes in the measured relaxation modulus can lead to arbitrarily large changes in the determined relaxation spectrum. In remedy, some reduction of the set of admissible solutions or appropriate regularization of the original problem can be used. Here, both techniques are used simultaneously.

Models
The modified spectrum is introduced: where the upper index of H M (v) means 'modified'. Then, (3) can be rewritten as i.e., the modulus is directly the Laplace integral of the spectrum H M (v).
is the space of real-valued squareintegrable functions on the interval (0, ∞). The respective sufficient conditions are given by Theorem 3 in [41]. Assume that the set of the linearly independent orthonormal functions {h 0 (v), h 1 (v), h 2 (v), . . .} form a basis of the space L 2 (0, ∞). Thus, the modified relaxation spectrum can be expressed as where the Fourier coefficients are [42] For practical reasons and in order to reduce the set of admissible solutions, it is convenient to replace the infinite summation in Equation (6) with a finite one of K terms, i.e., to approximate the relaxation spectrum H M (v) by a model of the form where the lower index of H M K (v) is the number K of model summands. Then, using (5), the respective model of the relaxation modulus is described by: where the functions Note, that function φ k (t) is Laplace transform of h k (v) for real argument t, i.e., with the notation L[ f (v)] used for Laplace transform. For basis functions, h k (v) applied in the developed algorithms, the components φ k (t) of the relaxation modulus model G K (t) are given by analytical or recursive formulas. This avoids quadrature errors occurring in the numerical calculation of the integrals (9). The following special functions are considered as basis functions: Laguerre, Legendre, Chebyshev and Hermite. The respective basis functions h k (v) and φ k (t) are described in Sections 3.3-3.6. All basis functions depend on one parameter-time-scaling factor α.

Identification Problem
Identification consists of selecting, within the given class of models defined by (7), (8) such a model, which ensures the best fit to the measurement results. Suppose a certain identification experiment (stress relaxation test [4,6,7]) performed on the specimen of the material under investigation resulted in a set of measurements of the relaxation modulus It is assumed that the number of measurements N ≥ K. As a measure of model (8) accuracy, the square index is taken where g K = g 0 · · · g K−1 T is a K-element vector of unknown coefficients of the models (7) and (8). Using the vector-matrix notation the identification index (10) can be rewritten in compact form as where · 2 denotes the square norm in the real Euclidean space R N . Thus, the optimal identification of relaxation spectrum in the class of functions defined by (7) and (8) consists of solving, with respect to the model parameter g K , the following least-squares problem: The matrix Φ N,K is usually ill-conditioned. Thus, the optimization problem (13) is still, like the original problem of solving Fredholm's equation of the 1st kind (3), incorrectly posed in the sense of Hadamard. Consequently, the solution of (13) is not unique, i.e., there exist many optimal model parameters minimizing the identification index Q N (g K ) (12). However, even the normal (with the lowest Euclidean norm) solution of (13) is noncontinuous and unbounded function of the measurement vector − G N . This means that when the data are noisy, even small changes in − G N would lead to arbitrarily large artefacts in any optimal model parameter. To deal with the ill posedness, the Tikhonov regularization method is used, as presented below.

Regularization
Regularization aims to replace the ill-posed problem with a nearby well-posed problem. Tikhonov regularization [43] strives to stabilize the computation of the least-squares solution by minimizing a modified square functional of the form: where λ > 0 is a regularization parameter. The above problem is well-posed; that is, the solution always exists, is unique, and continuously depends on both the matrix Φ N,K and on the measurement data − G N . The parameter vector minimizing (14) is given by: where I K,K is K dimensional identity matrix. The choice of regularization parameter λ is crucial to identify the best model parameters. Here, we apply the generalized cross-validation GCV [39,44], which does not depend on a priori knowledge about the noise variance. The GCV technique relies on choosing, as a regularization parameter, λ, which minimizes the GCV functional defined by [44] where the matrix are the residual vector for the regularized solution (15); tr[Ξ(λ)] denotes the trace of Ξ(λ). The problem of choosing the optimal regularization parameter has a unique solution and the resulting parameter differs the least from the normal solution of problem (14) that we would obtain for the ideal (not noise corrupted) measurements of the relaxation modulus [44].

Algebraic Background
Formula (15) is generally unsuitable for computational purposes. The singular value decomposition (SVD, [45]) technique will be used. Let SVD of the N × K dimensional matrix Φ N,K take the form [45]: where Σ = diag(σ 1 , . . . , σ r , 0, . . . , 0) R N,K is diagonal matrix containing the non-zero singular values σ 1 , . . . , σ r of the matrix Φ N,K [45], matrices V ∈ R K,K and U ∈ R N,N are orthogonal and r = rank(Φ N,K ) < N. Taking advantage of the diagonal structure of Σ and the matrices V and U orthogonality, it may be simply proved that the regularized optimal parameter − g λ K (15) is given by where K × N diagonal matrix Λ λ is as follows: Using SVD (18) and introducing N dimensional vector Y = U T − G N , the GCV function (16) can be expressed by a convenient analytical formula will consider several characteristics, as follows: (a) The model of the relaxation spectrum that we would obtain on the basis of ideal (undisturbed) measurements of the relaxation modulus: where ∼ g λ K is the vector of regularized solution of (14) for noise-free measurements of relaxation modulus G N = G(t 1 ) · · · G(t N ) T ; c.f., Equation (11) (b) The model of the relaxation spectrum that we obtain on the basis of the normal (13) for noise measurements of the relaxation modulus: where Φ † N,K is the Moore-Penrose pseudoinverse [48] of matrix Φ N,K , and  (13) for noise-free measurements of the relaxation modulus: Most of the results are formulated in terms of algebraic tools of the algorithm, i.e., SVD decomposition (18) of the matrix Φ N,K . Such an analysis enables a deeper insight into the properties of the algorithm and the resulting model. It shows not only the influence of the regularization parameter and measurement errors, but also the impact and significance of the selection of the basis functions, including their parameters and measurement points t i , on which the singular values σ i of the matrix Φ N,K depend.

Smoothness
The purpose of Tikhonov regularization relies on stabilization of the resulting vector where · 2 means the square norm, both in the real Euclidean space as well as in L 2 (0, ∞).
Therefore, the smoothness of the optimal solution − g λ K of discrete problem (14) guarantees that the fluctuations in the respective spectrum of relaxation, in particular the resulting spectrum of relaxation (22), are also bounded. In view of the above, due to orthonormality of the elements the basis system {h k (v)}, the function , optimal in the sense of the square identification index Q N (g K ) (10) of the bounded norm.
For any regularized   (20) and based on (28), the model smoothing efficiency can be evaluated by the following relation: which holds for an arbitrary regularization parameter λ > 0, where the spectrum − H N K (v) is given by (26). The last equality in (29) (26) and (28), the above result can be derived directly from the following inequality, proved in [28]: (29) illustrates the mechanism of stabilization. The following rule holds: the greater the regularization parameter λ is, the more highly bounded the fluctuations of the spectrum

Convergence
Relaxation spectrum (22) is only an approximation of that spectrum, which can be obtained in the class of models (7) by direct minimization (without regularization) of the quadratic index Q N (g K ) (12) for noise-free measurements, i.e., the approximation of the function T is vector of measurement noises, based on (28); the diagonal structure of Λ λ − Σ † and using the Schwarz inequality [49], the following bound of the relaxation spectrum approximation error can be derived: The inequality (30) yields that the accuracy of the spectrum approximation depends both on the measurement noises and the regularization parameter and on the singular values σ 1 , . . . , σ r of the matrix Φ N,K , which, in turn, depend on the selection of the basis orthogonal functions h k (v). Using (30), the regularized vector − g λ K converges to the noise-free normal solution ∼ g N K linearly with respect to the norm z N 2 , as λ → 0 and z N 2 → 0 , simultaneously. Therefore, the upper bound in (30) guarantees that the spec- at which they are both continuous, as λ → 0 and z N 2 → 0 , simultaneously.

Noise Robustness
The influence of disturbances in the measurements of the relaxation modulus on the regularized solution − g λ K was analyzed in detail in a two-part paper [27,50]. From Property 2 in [50], the following inequalities result: Thus, the regularized vector

Legendre Model
Let us assume a basis function with the time-scaling factor α, where P k (x) is Legendre polynomials [51][52][53] defined by Rodrigue's formula The polynomials P k (x) form a complete set of orthonormal basis in the interval [−1, 1] with the weight (2k + 1)/2 [51,53]. Thus, using the substitution x = 1 − 2e −2αv , it is easy to observe that the functions h k (v) defined by (31) form a complete orthonormal basis in L 2 (0, ∞) [49]. The relaxation modulus basis functions φ k (t) (9) are as follows: where the product ∏ p i=0 x i is equal to 1 when p < 0. The proof by induction is presented in Appendix A.1. The above formula can be equivalently expressed in recurrent form as Five first basis functions h k (v) are shown in Figure 1a,b for two different values of the time-scaling factor α. Figure 1c,d show the related φ k (t) functions. From the last figure, it is seen that the basis functions for the relaxation modulus model are in good agreement with the real relaxation modulus obtained in the experiment.

Laguerre Model
The Laguerre polynomials can be defined via Rodrigue's formula [51,54]: where > 0 is a time-scaling factor [55]. The continuous Laguerre function is the product of the Laguerre polynomial and the square root of the exponential weight function [56], i.e., The Laguerre functions form a complete orthonormal basis in (0, ∞) [51,56]. In Appendix A.2, the following formula is derived for the modulus basis functions: The above formula is given by Wang and Cluett [55], but as there are several definitions of the Laguerre functions in the literature, and, as a result, several formulas of the Laplace transforms, for example, in [57], the derivation of (35) is given in Appendix A.2 to avoid doubts.
A few first basis functions ℎ ( ) are shown in Figure 2a

Laguerre Model
The Laguerre polynomials can be defined via Rodrigue's formula [51,54]: where α > 0 is a time-scaling factor [55]. The continuous Laguerre function is the product of the Laguerre polynomial and the square root of the exponential weight function αe −αv [56], The Laguerre functions form a complete orthonormal basis in L 2 (0, ∞) [51,56]. In Appendix A.2, the following formula is derived for the modulus basis functions: The above formula is given by Wang and Cluett [55], but as there are several definitions of the Laguerre functions in the literature, and, as a result, several formulas of the Laplace transforms, for example, in [57], the derivation of (35) is given in Appendix A.2 to avoid doubts.
A few first basis functions h k (v) are shown in Figure 2a

Chebyshev Model
The Chebyshev polynomials of the first kind defined by the recursion relation [58,59]: starting with are orthogonal in the interval [−1, 1] with the weight function (1 − ) ⁄ [58]. Specifically, Thus, using the substitution = (1 − 2 ), it is easy to demonstrate that the set of functions with the first function defined as Basis functions Basis functions (34) and φ k (t) (35) of the Laguerre model for two time-scaling factors:

Chebyshev Model
The Chebyshev polynomials of the first kind defined by the recursion relation [58,59]: are orthogonal in the interval [−1, 1] with the weight function 1 − x 2 −1/2 [58]. Specifically, Thus, using the substitution x = 1 − 2e −2αv , it is easy to demonstrate that the set of functions with the first function defined as form an orthonormal basis in the space L 2 (0, ∞). Here, as previously, α is a positive timescaling factor. The relaxation modulus basis functions φ k (t) (9) are described by a useful recursive formula and for k = 0, 1, 2 are given by: where Γ(n) is the Euler's gamma function [60]. The proof is given in Appendix A.3, where two alternative formulas (A11) and (A12) for φ 1 (t) and φ 2 (t), respectively, are also derived. A few first basis functions h k (v) are shown in Figure 3a,b for two values of the factor α; the corresponding functions φ k (t) are plotted in Figure 3c,d. An earlier version of the model was presented in the paper [61]. form an orthonormal basis in the space (0, ∞). Here, as previously, is a positive timescaling factor. The relaxation modulus basis functions ( ) (9) are described by a useful recursive formula and for = 0,1,2 are given by: where ( ) is the Euler's gamma function [60].

Hermite Model
The Hermite functions defined as [52,62,63] where α > 0 is time-scaling factor and H k (x) are the Hermite polynomials, which satisfy recursion formula [52]: with the initial constitute an orthonormal system in the space L 2 (−∞, ∞) [52,63]. The relaxation modulus basis functions φ k (t) (9) are described by the recursive formula: and for k = 0, 1 are given by and where the complementary error function er f c(x) is defined by [64]: The derivation of the above formulas is given in Appendix A.4. The initial values H k (0) of Hermite polynomials are specified by [52]: A few first basis functions h k (v) are shown in Figure 4a,b for two factors α. The related functions φ k (t) are plotted in Figure 4c,d. The function φ 0 (t), and, in consequence of the recursive Formula (47), all functions φ k (t), depend on the exponential multiplier e t 2 /(2α 2 ) rapidly moving towards infinity. The following asymptotic properties were proved in Appendix A.5 for any α > 0: However, in numerical computations, the limited values of φ k (t) can be guaranteed only for t ≤ t upp , where t upp depends on the maximal real number accessible in the computing environment. For example, in Matlab, the largest finite floating-point number in IEEE double precision realmax = 2 − 2 −52 ·2 1023 ∼ = 1.7977·10 308 . Thus, in view of (48), the range of numerical applicability of the Hermite model in the time domain, determined by the inequality e t 2 /(2α 2 ) ≤ realmax, is as follows

Smoothness
Since the basis functions ℎ ( ) of the Hermite model presented above form an orthonormal basis of the space (−∞, ∞) of square integrable functions on (−∞, ∞), for Hermite model estimation (28) can be replaced by Therefore, for the Hermite model, the algorithm may (but does not have to) provide a stronger limitation of the fluctuation in the determined relaxation spectrum ( ) than for the other models.

Choice of the Basis Functions
In the models proposed above, the parameter > 0 is the time-scaling factor. The following rule holds: the lower the parameter is, the shorter the relaxation times are, i.e., the greater the relaxation frequencies are. The above is illustrated by Figures 1-4. Through the optimal choice of the scaling factor, the best fit of the model to the experimental data can be achieved. However, in practice, a simple rough rule for choosing the factor , based on the comparison of a few first functions from the sequence

Smoothness
Since the basis functions h k (v) of the Hermite model presented above form an orthonormal basis of the space L 2 (−∞, ∞) of square integrable functions on (−∞, ∞), for Hermite model estimation (28) can be replaced by Therefore, for the Hermite model, the algorithm may (but does not have to) provide a stronger limitation of the fluctuation in the determined relaxation spectrum − H M K (v) than for the other models.

Choice of the Basis Functions
In the models proposed above, the parameter α > 0 is the time-scaling factor. The following rule holds: the lower the parameter α is, the shorter the relaxation times are, i.e., the greater the relaxation frequencies are. The above is illustrated by Figures 1-4. Through the optimal choice of the scaling factor, the best fit of the model to the experimental data can be achieved. However, in practice, a simple rough rule for choosing the factor α, based on the comparison of a few first functions from the sequence {φ k (t)} for different values of α with the experimentally obtained function − G(t i ), is quite enough. In the same manner, the number K of the series G K (t) (8) elements can be initially evaluated. This rough selection strategy of the model parameters was used in the examples presented below. Thus, the choice of the number K and the parameter α must be carried out a posteriori, after the preliminary experiment data analysis.
The ranges of applicability of the four classes of models described above in the relaxation times t domain and the relaxation frequencies v domain for different values of α are summarized in Table A1 in Appendix A.6. It was assumed that the range of applicability for times is determined by the value of time t, for which the first K = 11 basis functions φ k (t) no longer permanently exceed, i.e., for any θ > t, ε = 0.5% of its maximum value. Specifically, where Similarly, the range of applicability for the relaxation frequencies was defined on the basis of the variability in the basis functions h k (v). Here, with h kmax defined by In view of the problems described above concerning the numerical determination of the basis functions only for t ≤ t upp , with t upp defined in (54), for the Hermite model, ε = 0.0212 was assumed.

Example 1
Consider viscoelastic material of relaxation spectrum described by Gauss-like distribution [11] H The corresponding modified spectrum H M (v) (4) is: and, therefore, using (5), the 'real' relaxation modulus is where er f c is defined by (50). In the experiment, N = 1000 sampling instants t i were generated with the constant period in the time interval T = [0, 32] seconds selected in view of the course of the modulus G(t) (59). Additive measurement noises z(t i ) were selected independently by random choice with uniform distribution on the interval [−0.005, 0.005]Pa, i.e., maximally 6.1% of the mean value of G(t) in the interval T defined as the average value of the integral of G(t) over T , which is equal to 0.0820Pa. The time-scaling factors α are selected by comparison for different values of α a few first functions φ k (t) with the experiment results − G(t i ). Only for the Chebyshev model, the rough selection of α required several attempts; for the remaining classes of models, it was enough to review the data from Table A1. The basis functions h k (v) and φ k (t) were simulated in Matlab R2022a using special functions erfc, legendreP, chebyshevT, and hermiteH. For the singular value decomposition procedure, svd was applied. For K = 6, 8, 9, 10,11, and sometimes also for K = 12, the regularization parameters λ GCV were determined and are given in Table 1.
Next, the vectors of optimal model parameters − g λ K = − g λ GCV K (19) were computed and are given in Table A2 in Appendix B. Only for the Laguerre models, some elements of the vectors of optimal parameters are negative; for the remaining classes of models, all the parameter vectors − g λ GCV K are positive. The optimal modified spectra of relaxation frequencies − H M K (v) (22) and the 'real' spectrum H M (v) (58) are plotted in Figure 5 for selected values of K, while in Figure 6, the spectra − H K (v) (23) and H(v) (57) are presented. The optimal models of the relaxation modulus G K (t) computed for g K = − g λ GCV K according to (8) are plotted in Figure A1 in Appendix B, where the measurements − G(t i ) are also marked. Since the optimal models G K (t) (8) and is given in the last column of Table 1.

Example 2
Consider the spectrum of relaxation times introduced by Baumgaertel, Schausberger and Winter [36,37], which is known to be effective in describing polydisperse polymer melts [17,34], with the parameters [34]: β 1 = 6.276·10 4 Pa, β 2 = 1.27·10 5 Pa, τ c = 2.481s, τ max = 2.564·10 4 s, ρ 1 = 0.25 and ρ 2 = −0.5. The related spectra of relaxation frequencies H(v) = H 1 v and H M (v) (4) are well posed for v > 0. The modified spectrum is described by and depicted in Figure 7; the corresponding 'real' relaxation modulus G(t) is defined by (5). In the experiment, N = 1000 time instants t i were sampled according to the square rule t i = ∆t(i − 1) 2 + 60 s, with parameter ∆t = T/(N − 1) 2 , where T = 10 7 s, in the time interval T = [0, T] is selected in view of the course of the modulus G(t). Due to the numerical problems described above related to determining the basis functions φ k (t), Equations (47)- (49), for the Hermite model the experiment was simulated in a shorter time interval T = 10 6 s using the same sampling formula. Additive measurement noises z(t i ) were selected independently by random choice with uniform distribution on the interval [−0.0005, 0.0005]MPa. Here, for the selection of parameter α, which guarantees a satisfactory accuracy of the modulus approximation, several or even a dozen or so attempts were necessary. This means that it will be advisable to further extend the algorithm by the level of optimal selection of the time-scaling factors. For K = 6, 8, 10,12,14, and also for K = 16, the regularization parameters λ GCV were determined and are given in Table 2.
Only for the Laguerre and Legendre models, the accuracy of the modulus approximation , Equation (10), is comparable to that obtained in Example 1.
Thus, the vectors of optimal model parameters 19) for Legendre and Laguerre models are given in Table A3 in Appendix B; for the remaining models, they were omitted.
The optimal spectra − H M K (v) (22) and the 'real' spectrum H M (v) (58) are plotted in Figure 7 for Laguerre and Legendre models; a logarithmic scale is used for the frequencies v and a linear scale for the spectrum. For Legendre and Laguerre models, the optimal models G K (t) computed for g K =   (22) was estimated by the integral square errors ERR, defined in (60) and given in Table 2, only for Legendre and Laguerre models. For the Chebyshev and Hermite models, a satisfactory quality of approximation was not obtained, despite the research carried out for a wide range of α values. The best found α and related identification indices Q N − g λ GCV K are given in Table 2, but other indices are omitted.
The optimal Legendre and Laguerre models G K (t) (8) were well fitted to the experimental data, see Figure 8a Table 2). In particular, the course of the G K (t) model from Figure 8c shows that it is not too few components of the model but the properties of the Chebyshev model that make it ineffective for approximating the relaxation spectrum H M (v) (61).

Example 2
Consider the spectrum of relaxation times introduced by Baumgaertel, Schausberger and Winter [36,37], which is known to be effective in describing polydisperse polymer melts [17,34], with the parameters and depicted in Figure 7; the corresponding 'real' relaxation modulus ( ) is defined by (5). In the experiment, = 1000 time instants were sampled according to the square rule = ∆ ( − 1) + 60 s, with parameter ∆ = ( − 1) ⁄ , where = 10 , in the time interval = [0, ] is selected in view of the course of the modulus ( ). Due to the numerical problems described above related to determining the basis functions ( ), Equations (47)- (49), for the Hermite model the experiment was simulated in a shorter time interval = 10 using the same sampling formula. Additive measurement noises ( ) were selected independently by random choice with uniform distribution on the interval [−0.0005, 0.0005] . Here, for the selection of parameter , which guarantees a satisfactory accuracy of the modulus approximation, several or even a dozen or so attempts were necessary. This means that it will be advisable to further extend the algorithm by the level of optimal selection of the time-scaling factors. For = 6,8,10,12,14, and also for = 16, the regularization parameters were determined and are given in Table 2. Only for the Laguerre and Legendre models, the accuracy of the modulus approximation measured by , Equation (10), is comparable to that obtained in Example 1.
Thus, the vectors of optimal model parameters = (19) for Legendre and Laguerre models are given in Table A3 in Appendix B; for the remaining models, they were omitted. The optimal spectra ( ) (22) and the 'real' spectrum ( ) (58) are plotted in Figure 7 for Laguerre and Legendre models; a logarithmic scale is used for the frequencies and a linear scale for the spectrum. For Legendre and Laguerre models, the optimal models ( ) computed for = are plotted in Figure 8, where the measurements ̅ ( ) are also marked and logarithmic scale is applied for time axis. The norms , ‖ ( )‖ and the identification index are given in Table 2. For the    , Equation (10), the errors ERR (60) of the relaxation spectrum models.

Remarks
The example orthogonal basis functions for relaxation spectrum models have been assumed as the products of exponentials and Legendre, Laguerre, Chebyshev and Hermite polynomials. For Legendre and Laguerre models, the basis functions φ k (t) of the relaxation modulus are rational functions; for the Chebyshev model, they are determined by the quotient of the Euler's gamma functions, which is a generalization of the factorial of a non-negative integer for the no-integer argument, while, for the Hermite model, these functions are based on the complementary error function. From Figures 1c,d and 2c,d, it is evident that the basis functions φ k (t) for the relaxation modulus of Legendre and Laguerre models are in good agreement with the real relaxation modulus obtained in the experiment. However, Figures 3c,d and 4c,d show that, for the Chebyshev and Hermite models, so good agreement is not achieved for individual functions φ k (t). Hence, a significantly worse fit of the Chebyshev and Hermite models to the measurement data for the spectrum from Example 2, which has a much wider range of relaxation frequencies, results.
Both examples show that, generally, increasing the number of model summands improves the model quality, provided that the assumed series can provide a good approximation of the relaxation modulus. If a given series of special functions is not suitable for approximation of the relaxation modulus for a given material, then, as shown by the research conducted for the Chebyshev and Hermite models in the second example, increasing the number of series components does not improve the quality of the model.

Remarks
The example orthogonal basis functions for relaxation spectrum models have been assumed as the products of exponentials and Legendre, Laguerre, Chebyshev and Hermite polynomials. For Legendre and Laguerre models, the basis functions ( ) of the Both examples show that, generally, increasing the number of model summands improves the model quality, provided that the assumed series can provide a good approximation of the relaxation modulus. If a given series of special functions is not suitable for approximation of the relaxation modulus for a given material, then, as shown by the research conducted for the Chebyshev and Hermite models in the second example, increasing the number of series components does not improve the quality of the model.

Conclusions
In this paper, a class of algorithms for the relaxation spectrum identification, which combines the technique of an expansion of a function into a series in an orthonormal basis with the least-squares regularized identification, was derived. It was demonstrated that due to the choice of an appropriate special functions (Legendre, Laguerre, Chebyshev,

Conclusions
In this paper, a class of algorithms for the relaxation spectrum identification, which combines the technique of an expansion of a function into a series in an orthonormal basis with the least-squares regularized identification, was derived. It was demonstrated that due to the choice of an appropriate special functions (Legendre, Laguerre, Chebyshev, Hermite) as the basis functions for the relaxation spectrum model, the basis functions of the relaxation modulus model are given by compact analytical or recursive formulas. Due to the choice of orthonormal basis, the smoothness of the vector of the optimal model parameters implies equivalent smoothness of the fluctuations in the model of the relaxation spectrum.
The proposed approach based on the expansion of the relaxation spectrum into a function series can be applied for arbitrary basis functions. The proven convergence and noise robustness properties of the optimal models will be retained, but the smoothing of the spectrum model will require separate analysis. As part of further research, the algorithm can be extended with a superior level of optimal selection of the time-scaling factor, so as to obtain a better fit to the measurement data. The presented scheme of the relaxation spectrum identification can be easily modified for the retardation spectrum recovery from the creep compliance measurements obtained in the standard creep test. Mathematical induction will be used. In the proof, the following recurrence representation of the Legendre polynomials [52,53]: is used, where the two first polynomials are as follows and In the base case of the proof by induction, for k = 0 and k = 1, Formula (32) follows immediately from definition (9) and properties (A2), (A3). For k = 0, by (9) (31) and (A2) the basis function φ k (t): , while for k = 1, applying (A3), we obtain: .
Appendix A.2 Derivation of Formula (35) To derive the basis function φ k (t) (35), we start with the alternative formula of the Laguerre polynomials L k (v) (33) given by [65]: An easy consequence of (A6) is where, using the well-known formula for the Laplace transform the basis functions φ k (t) defined by (9) are given by for k = 0, 1, 2, . . .. Hence, by binomial theorem, we have where Formula (35) follows. This formula is known in [55] and coincides with that given in [66].
Comparing the first integral in Formula (A13) and the integral φ 1 (t), it is easy to notice that Since by (50), we have der f c t/ √ 2α Equations (A14) and (48) yield where the next formula follows and, in view of (48), Formula (49) results. The proof of the recursive Equation (47) is based on (45), which, by definitions (44) and (9)  and taking into account the known property of Hermite polynomials [52]: and definitional Formulas (44), (9), can be expressed as   (55) and v app (56) of the respective intervals 0, t app and 0, v app are given.     are expressed in MPa·s 1/2 .

Model
Legendre model